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We study the scattering of light (null geodesies) by two fixed extreme Reissner-Nordstrom 
black holes, in which the gravitational attraction of their masses is exactly balanced with the 
electrostatic repulsion of their charges, allowing a static spacetime. We identify the set of 
unstable periodic orbits that form part of the fractal repeller that completely describes the 
chaotic escape dynamics of photons. In the framework of periodic orbit theory, the analysis 
of the linear stability of the unstable periodic orbits is used to obtain the main quantities 
of chaos that characterize the escape dynamics of the photons scattered by the black holes. 
In particular, the escape rate which is compared with the result obtained from numerical 
simulations that consider statistical ensembles of photons. We also analyze the dynamics 
of photons in the proximity of a perturbed black hole and give an analytic estimate of the 
escape rate in this system. 

PACS numbers: 05.45.-a, 04.70.Bw 



I. INTRODUCTION 



In general relativity, the nonexistence of an absolute time parameter introduces new aspects 
in the characterization of chaos with respect to the well-known Newtonian dynamics Some 
relativistic systems in which the existence of chaos has been reported include charged particles 
in a magnetic field interacting with gravitational waves, spinning particles orbiting rotating and 
nonrotating black holes, gravitational waves from spinning compact binaries, as well as particles 
in Majumdar-Papapetrou geometries. 

Most of the studies of chaos around black holes have focused on the analysis of the qualitative 
changes in the dynamics of an isolated black-hole spacetime caused by a small perturbation due to 



external mass distributions 110,018, ^vitationai waves fl, spin-orbit and spin-spin coupling 
a], or magnetic fields [9[]. Recently, the projects to startup operative ground-base gravitational 
wave detectors (LIGO, VIRGO, GEO600, TAMA300, AIGO) and a planned laser interferometer 



space antenna (LISA) [lOj], which will be able to detect the signals from gravitational wave sources, 
such as inspiralling compact binary systems of neutron stars or black holes, have increased the 
interest in the presence of chaos in the dynamics of binary black holes and its effects on the 



outgoing gravitational radiation 



111. 



In this work we use the Majumdar-Papapetrou metric [12| to analyze the dynamics of photons 
in the gravitational field of two extreme Reissner-Nordstrom black holes that are fixed in space 
due to the balance of their gravitational attraction and electrostatic repulsion. Although it is 
unlikely that this metric describes any astrophysical system, since in real universe black holes tend 
to rotate around their center of mass producing gravitational waves and do not possess overall 
electric charge, the chaotic scattering of photons in the Majumdar-Papapetrou static spacetime of 
nonrotating black holes with extreme electric charge still provides an interesting formal model that 
can be used to illustrate many of the features expected in more realistic systems. 

In [2| Chandrasekhar studied the scattering of radiation by two extreme Reissner-Nordstrom 
black-holes that are at finite distance apart. In contrast to the two center problem in Newtonian 
gravitation, in general relativity this two center gravitational problem is generally not integrable. 
An appendix in displays a set of null geodesies (photon's trajectories) in the meridian plane of 
the system and concludes that the dynamics is probably not separable. 

Shortly after, Contopoulos [3J systematically studied the set of photon trajectories with zero 
angular momentum component along the axis that goes through the black-holes. In that case the 
motion is confined to a plane (the meridian plane). The study concluded that there are three types 
of non periodic motions. There are orbits that fall into one of the black-holes, with mass M\ or 
mass M2, these are orbits of type / and II, and there are orbits that escape to infinity, these are 
orbits of type III. The orbits of different types are separated by orbits that tend asymptotically 
to three main types of unstable periodic orbits. One kind of periodic orbits goes around one of the 
black-holes, either Mi or M<i- A third type of unstable periodic orbits goes around the two black- 
holes. Between two non periodic orbits of two different types there are orbits of a third type. In a 
related work Contopoulos et at. [13] analyze in detail the different types of periodic orbits in the 
two fixed center system and compare families of these orbits in the Newtonian and the relativistic 
problems. In this system chaos appears explicitly as the initial conditions of these types of orbits 
form a Cantor set. Further numerical evidence indicated that all the photon periodic orbits are 



3 



unstable. In the same direction it was investigated the phase space trajectories in a multi-black-hole 
spacetime 4|, it was found that the chaotic geodesies are well described by Lyapunov exponents. 
All these works strongly indicate that the scattering of photons by two Reissner-Nordstrom black- 
holes held fixed is chaotic and that the set of periodic orbits is unstable. The geometric analysis 
of the flow also showed the chaotic behavior of the relativistic null-geodesic motion in the two 
black hole spacetime [141 ] . More recently Contopoulos et al. have studied in detail the asymptotic 
curves from the periodic orbits, their homoclinic and heteroclinic intersections and the basins of 



attraction of two black holes 



Our aim is to give a full description of the chaotic escape dynamics of photons in the two 
black-hole system in terms of the different physical indicators of chaos. We shall analyze the linear 
stability of the dynamics, in particular the unstable periodic orbits, and from that information 
we evaluate the escape rate associated with the equations of motion using the well known trace 



formula for hyperbolic flows of Cvitanovic and Eckhardt [10]. It turns out that the escape rate 
is given by the leading eigenvalue of the spectrum of the evolution operator [17]. In addition, we 
shall compare our results with numerical simulations that consider the time evolution of statistical 
ensembles of photons. 

In Section II we describe the Chandrasekhas's model to analyze the photon dynamics under the 
gravitational field of two extreme Reissner-Nordstrom black holes. In Section III we discuss the 
linear stability of the dynamics and the general methods that we will consider to obtain the escape 
dynamics of photons. Section IV introduces the formalism to study the time evolution of statistical 
ensembles of photons and their escape dynamics. In Section V we present the set of periodic orbits 
of the system in the meridian plane {L z = 0) and give their periods and stretching factors. In 
this section we also illustrate the fractal repeller associated with the dynamics of photons in the 
two black-hole field. In Section VI we calculate the main quantities of chaos derived from the 
analysis of the linear stability of the unstable periodic orbits. The escape rate obtained from the 
periodic orbit theory is compared with the value obtained from numerical simulations that consider 
statistical ensembles of photons. In Section VII we analyze the dynamics in the proximity of a 
perturbed black hole and give an analytic estimate of the escape rate in this system. In Section 
VIII the main conclusions are put together. 
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II. DYNAMICAL EQUATIONS 

The model studied by Chandrasekhar in Q] derives from the solutions of the Einstein-Maxwell 
equations to describe a problem analogous to the Newtonian arrangement of charged mass-points 
in which the mutual Coulomb repulsions are balanced with the gravitational attraction. These 
solutions are known as the Majumdar-Papapetrou solutions [ijj]. They are obtained from an static 



solution of the Einstein-Maxwell equations 



la ]. The metric of the Majumdar-Papapetrou solution 



of the Einstein-Maxwell equations is given by 

ds 2 = dt 2 /U 2 -U 2 (dx 2 + dy 2 + dz 2 ), (1) 

where (x, y, z) are the spatial coordinates. The function U, which depends only on the spatial 
coordinates, is a solution of the three-dimensional Laplace's equation 

^ ^ (^dx 2 dy 2 dz 2 ^ ^ ^ 
Hartle and Hawking [19] showed that for a function U of the form 

N M- 

C/ = l + E^ (3) 

i=i r < 

being r» = y/(x — xi)' 2 + (y — yi) 2 + (z — Zi) 2 , the Majumdar-Papapetrou solution corresponds to 
a system of N extreme Reissner-Nordstrom black holes with horizons at (xi,yi, Zi) and with masses 
equal to their charges, Mi = Qi > 0. The metric is everywhere regular except at the black-hole 



locations where there are coordinate singularities, as argued in 

We will consider a configuration of two black holes located at (0, 0, ±Zbh)- Using geometrized 
units (the speed of light in vacuum c = 1 and the gravitational constant G = 1) the function U 
takes the form 

n Mi M2 

(x 2 + y 2 + (z- z bh ) 2 y/ 2 + (x 2 + y 2 + (z + z bh ) 2 y/ 2 ' 1 ' 

The Lagrangian C associated with the metric ([TJ is defined by 

C = ^ r2 -^(x 2 + y 2 + z 2 ), (5) 

where the dot denotes the derivative with respect to the affine parameter [i. The Hamiltonian 7i 
can be expressed as 

n= l -u 2 p 2 - l -u- 2 (p 2 x + p 2 y +pl), (6) 
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being the Hamilton's equations that dictate the geodesic motion associated with the metric ([T]): 

^ Q tt-2 frr\ 

a = - = -U* Pa (7) 
P« = ^ = ~\d«{U 2 - U- 2 P 2 ) (a = x,y,z). 

According to the Hamilton's equations pt = £ = constant and angular momentum along the 
z-axis L z = xp y — yp x = constant. In the case of photons it follows also that TL = U 2 £ 2 /2 — 
U~ 2 (pl +P 2 +p1)/2 = 0, [J, and therefore p 2 + p 2 + p 2 = P 2 = U 4 £ 2 . Hence the problem scales 
with £ and without loss of generality we can consider £ = 1. 

The equations of motion with respect to the killing time t can be written as: 

, da .du 4 

a = Tt =a Tt = - u Pa (8) 

P ' a = ~a¥ = Pa ~a% = ~ 2U ~ ldaU (a = x,y,z). 

The sets of equations (|7|) and ([8]) are related by the transformation 

± = U 2 ^- (9) 
du dt w 

between the affine parameter u and the killing time t. Although the solutions of these two systems 
induce different flows, and $ respectively, and therefore lead to different values in some prop- 
erties that depend on the flow, the spatial distribution of their corresponding set of orbits do not 
change, and neither are modified the critical elements of the flow, in particular the periodic orbits. 

Due to the conservation of L z , the motion for L z = is restricted to a plane (the meridian 
plane). We shall take this simplification and study the dynamics in the meridian plane. 



III. LINEAR STABILITY 



We are particularly interested in the stability analysis of the solutions of the equations of motion, 
i.e. how is they behave after being slightly perturbed. To fix the notation and make the paper 
self contained, in this section we introduce the basic concepts of linear stability. First we shall 
discuss this problem in a general context. Then the application to the two black hole system is 
straightforward. 

Let us consider a Hamiltonian system with / degrees of freedom. The 2/-dimensional phase- 
space can be denoted by X = (X\, X2, ■ • • ,X 2 f), being (Xi, ■ ■ ■ , Xf) the generalized coordinates 
and (Xf+i, ■ ■ ■ , X 2 f) the conjugated momenta. The Hamilton's equations for a Hamiltonian TL can 
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be expressed as 



being X = dXj dr, with r the afhne parameter fi or the killing time t, and S a 2/ x 2/ antisymmetric 
matrix of the symplectic form 



£ = 



( ^ 

V" 10 / 



(ii) 



with and 1 the / x / null and unit matrices respectively. 

A solution X T of these differential equations is a curve in phase-space corresponding to a given 
initial condition Xq. The stability analysis of this solution involves the study of the time evolution 
of a small perturbation 5X with respect to X T . Considering Y T = X T + 5X, it follows that SX 
satisfies the initial value problem 

SX = L(r) • SX, (12) 

with L(r) =S-0| Xr . This is a linear initial value problem with time-dependent coefficients. For 
a given initial perturbation SXq it has the solution 

SX T = M(X ,t) -SX . (13) 

The If x 2/ matrix M(Xq,t) is usually referred as the fundamental matrix and gives the time 
evolution of an initial displacement SXq. The fundamental matrix satisfies the differential equation 

M(X ,t) = L(t)-M(X ,t), with M(X ,Q) = 1. (14) 

The behavior of the perturbation SX can be analyzed in terms of the Lyapunov exponents 
of X T , which measure the rate of exponential separation or approach of initially infinitely close 
trajectories. The Lyapunov exponent associated with the unit vector ej = SXj/\SXj\ along the 
direction of SXj is given by 

ApT ,e 7 -)= lim -ln|M(X ,r) -e 7 |. (15) 

J T— >0O f 

The number of Lyapunov exponents is equal to the phase-space dimension. Positive Lyapunov 
exponents indicate dynamical instability of trajectories in phase space, and therefore an extreme 
sensitivity to the initial conditions. 

For a periodic orbit X T , the Lyapunov exponents are directly related to the stability of this 
solution. The Lyapunov exponents can be degenerate with multiplicity mj, in the sense that 
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several of them have the same value. They are ordered as Ai > A2 > • • • X v with ^ rrn = 2f. 

i=l 

The fundamental matrix has a symplectic structure and therefore satisfies the relation M r SM = S 
[20I ]. Due to this property all the Lyapunov exponents in Hamiltonian systems are grouped in 
pairs of equal absolute value and opposite sign {Aj, — Xi}i=i t ... j. Furthermore, perturbations along 
the direction of the flow give a zero Lyapunov exponent, and due to the pairing property another 



21 



22]. Thus at least two 



one is also zero, which is a consequence of the conservation of energy 
exponents are zero. 

The matrix M(Xq, t) corresponding to a periodic orbit is also periodic with the same period. In 
general X T = X t+ t, with T = rT p (r=l,2,- • •) and T p the primitive period of the orbit. The linear 
stability of a periodic orbit is determined by the eigenvalues of M(Ao, r) over one primitive period, 
i.e. M(Xq, T p ). These eigenvalues, called stretching factors or stability eigenvalues, are obtained 
from the secular determinant 

det[M(X ,T p ) - Al] = 0. (16) 

Due to the symplectic property of the fundamental matrix the eigenvalues Aj are also grouped 
in pairs. From the previous discussion two of them are equal to one. One of them corresponds 
to perturbations along the periodic orbit and the other to perturbations that are perpendicular to 
the energy shell. The Lyapunov exponents associated with a periodic orbit are given in terms of 
the stretching factors as Aj = In |Aj|/T p . 

The stability of a particular motion is determined by the location of the stretching factor in the 
complex plane: |Aj| > 1 is a signature of instability and corresponds to the unstable manifold of 
the periodic orbit, |Aj| < 1 is related to the stable manifold and |Aj| = 1 corresponds to the center 
manifold of the orbit. 

A straightforward application, which is directly related to our problem, is the study of the 
dynamics of two-dimensional systems. In these systems a periodic orbit is called hyperbolic when 
A > 1 > A -1 > 0, hyperbolic with reflection when A < —1 < A -1 < 0, and elliptic when |A| = 1 



20|. The orbit stability changes whenever A crosses the unit circle, and generally a bifurcation 
may occur in which some periodic orbits can disappear or emerge. In autonomous systems the 
energy is the main bifurcation parameter. However, in the two black-hole system no bifurcation is 
expected as the photon dynamics scales with the energy. 

A periodic orbit is characterized by its primitive period and its stretching factors Aj. We shall 
see below how this information can be used to analyze in detail the escape dynamics of photons 
scattered by two black-holes. In general the Lyapunov exponents and the periods of the periodic 
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orbits are not invariant under time-transformations, but the stretching factor are invariant, as 



noted in 



, 6, Q, Qj- In particular Motter 



l24j |. In particular Motter [25[] showed that if we consider a transformation of 
the form dp = f(X)dt, the Lyapunov exponents transform according to Af = A'/ < / >t, where 
< / >t= limt-^»ao(l / 1) Jo f{X{t?))dlf. For a given periodic orbit p the Lyapunov exponents over 
a period in the /i-time and over a period T in the i-time are related by xj" =)£*/</>* 
with < / >* T = (l/r t )/ Jt f(X p (t'))dt' = Tp/Tt. From these relat ions it can be concluded that 

T T T T 

Aj M T M = X i t T which implies the invariance of the stretching factors, = A^*. 

IV. EVOLUTION OF THE DENSITY OF PHOTONS AND THEIR DECAY DYNAMICS 

The motion of a single photon in the meridian plane is chaotic and it is useful, in order to 
characterize its escape dynamics from the black holes, to consider statistical ensembles that provide 
a probabilistic description of the process. Therefore, instead of analyzing a single photon orbit, we 
consider an ensemble described by some density p. For a given Hamiltonian H the time evolution 
of a density p is dictated by the Liouville equation 

d T p = Lp, (17) 

with the Liouville operator L given by the Poisson bracket as L = {7i,p}. Assuming a time- 
independent Hamiltonian, the solution of this linear equation can be expressed as 

Pr = e% , (18) 

being po the initial density. From (|18p . the average over the ensemble of any observable A at time 
r can be computed as; 

(A) r = J dX A($ T X)po(X), (19) 

being $ T the flow that maps an initial condition Xq into X T , according to X T = <& t Xq. This last 
expression can also be written in the form 

(A) T = J dXdYA(X)5(X - $TY)po(Y), (20) 

where the Dirac distribution 5{X — & T Y) represents the conditional probability density of a tra- 
jectory being at position X at time r given the initial position Y. Indeed the Dirac distribution 
defines the kernel of the evolution operator. The equation (|20p defines two self adjoint operators 
as 

(A) T = (P^A\ P0 ) = (A\P T po) . (21) 
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The evolution of a probability density is ruled by the Frobenius- Perron equation 2l|, |22 



Pr{X) = P T p (X) = J dY5(X - $ T Y)p (Y) 



(22) 



and the time evolution of the observable is dictated by the Koopman operator 26| defined by 



A T (Y) = P^ T A{Y) = J dX5(X - $ T Y)A(X). 



(23) 



In the case of an invertible and conservative flow <3? T the Frobenius-Perron operator reduces to 
p T {X) = P T po(X) = po(&~ T X). In open systems the trajectories that are initially confined in a 
bounded region in phase space tend to escape towards infinity with an exponential escape rate. The 
number of particles iV(r) that remain inside the initial domain at instant r decays exponentially 
with time as N(t) ~ N(0)e~ JT , being 7 the so called escape rate. 

The leading eigenvalue of P T dominates the decay and determines the escape rate. The rest of 
the spectrum (resonances) describe further details of the dynamics and give other important time 
scales besides the one associated with the escape rate. Thus the spectrum of the Frobenius-Perron 
operator P T provides a way to describe the evolution of a set of trajectories characterized by a 
density p T (X). 

~ n □ 

The spectral theory of P T developed from the works of Koopman [26| and von Neumann [271 ]. 
It assumes that the evolution operator acts on a functional space of square integrable densities C? . 
In that case the evolution is unitary and the eigenvalues of P T are located on the unit circle. The 
spectrum of a chaotic system presents continuous components on the unit circle which describe 
the correlation function decay for these systems. New methods have been developed to obtain the 
eigenvalues of P T outside the unit circle and to analyze the escape process or relaxation dynamics 



28, 



291 ]. These methods, which are valid for systems where all the periodic orbits are unstable, 



have been successfu 
applications 



21 



22 



ly used in chaotic systems and the study of transport phenomena among other 
, 32]]. We shall apply them to the study of the escape dynamics of 



30 



31 



photons from the two black holes. 

One of the methods developed to obtain the spectrum computes the trace of the Frobenius- 



Perron operator. Here we shal 
literature on the subject 



21 



22 



outline this procedure and refer the reader to the specialized 



28, 



29 



331 ] for further details. 



The trace of P T can be formally written as 



TrP T 



dX5(X - $ T X). 



(24) 
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The contributions to the trace are due to the fixed points of the flow, which are given by the 
condition 

X = $ T X, (25) 

at some r. In general this equation presents two types of solutions: stationary points that satisfy 
(|25p for all r, and periodic orbits which satisfy the fix-point condition at a discrete set of values, 
given by r = rT p , with r = 1,2, We assume that the solutions of ([25]) are isolated, and 
therefore each stationary point or periodic orbit is locally unique. A periodic orbit may belong 
to a continuous family provide that there exists some continuous symmetry or some constant of 
motion. When this occurs the integration domain in (|24[) must be reduced to a less dimensional 
space until the periodic orbit is completely isolated. In time-independent Hamiltonian systems 
the periodic orbits are rarely isolated, indeed they tend to form continuous families as the energy 
changes. Hence the phase-space must be reduced considering the energy conservation 7i = E. 
Then the trace ()24p is formally given by 

Tr E P T = f d 2f - l x8{X-4> T E x), (26) 

JE 



when* <>y denotes the flow on the energy shell [nl If any additional symmetry is present in 
the system then further reductions are required. In the two black system we are considering the 
dynamics on the energy shell and we have removed the axial symmetry by choosing a particular 
meridian plane. Therefore all the periodic orbits of the flow can be treated as isolated. 

The tra ce (1261) can be written in terms of the unstable periodic orbits of the flow and their 



repetitions [la . |34|. In order to do so, the integral in (j26|) is resolved considering at each point a 
coordinate system with one axis fixed along the periodic orbit. The integral over the coordinate 
along the orbit is trivial, being related to the period of the orbit. The integration over the transverse 



coordinates takes into account the stability of the periodic orbit. The result for the trace is 



where the sums are over the primitive periodic orbits (p.p.o.) and their r repetitions. The matrix 
m p is derived from the fundamental matrix M(Xq,T p ) once the perturbations along the orbit and 
the perturbations perpendicular to the energy shell have been removed. Hence in a two-dimensional 
flow the matrix m p only contains the stretching factors whose modulus is different from one. These 
are precisely the ones that are related to the stable and unstable manifolds of the isolated unstable 
periodic orbit. 
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Formally the Laplace transform of the Frobenius-Perron operator gives the resolvent of the 
Liouville operator defined in (|17p . which using (I26D can be written as 

1 r°° d 

Tr ^= dTe- aT TrP = —In Z(a), (28) 

a - L Jo oa 



where Z{a) is the so called Selberg-Smale Zeta function 



given by 



*M = «p(- E E; | det(1 - m ai )- < 29 > 

\ p=p.p.o. r=l 1 v P' 1 / 
Denoting by A p and A" 1 the eigenvalues of m p , the zeta function can be expressed as 



_ aT \ fe+l 
e al p 



m= n n i-y ^ 

p=p.p.o k =0 \ \ ly P\ lv pJ 

where A p > 1 since all the periodic orbits are unstable in the two black-hole system. The spectrum 
of the system can be determined from the zeroes of the function Z(a), which is usually expressed 
as a product of the inverse Ruelle ^-functions, 

^) = c ~ V)cr VR 2 - V) • • • (3i) 

with the definition 

^-^v-S)- (32) 

The relaxation times 7$ associated with the dynamics are given by the poles of the resolvent of 
the Liouvillian, as Reai = rf 1 , or equivalently, by the zeroes of the zeta function or the zeroes 
of the inverse of the Ruelle (^-functions. The leading part of the spectrum is given by the poles of 
the first Ruelle £- function Co(c) (or zeroes of Co" 1 )- So we are mainly interested in the part of the 
spectrum that controls the long-time behavior of the dynamics. 

The set of unstable periodic orbits can be used to describe in fine detail the escape dynamics 
of photons from the black holes. The number of periodic orbits in the fractal repeller grows 
exponentially with the period, with a rate of proliferation that is given by the topological entropy 
htop- The amount of chaos in the system is given by the Kolmogorov-Sinai entropy per unit time 
h,KS, which measures the minimal data accumulation rate required to reconstruct a trajectory on 
the repeller without ambiguity. The repeller has a mean Lyapunov exponent A. In a Poincare 
section that is transverse to the flow on the repeller, the fractal repeller defines another fractal 
which is characterized by the partial generalized fractal dimensions. Two representative partial 
dimensions are the Hausdorff dimension dn, which characterizes the bulkiness of the repeller in 
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phase space, and partial information dimension dj that can be obtained from the Kolmogorov-Sinai 
entropy and the mean Lyapunov exponent as di = hxs/^ 3]- 

The different quantities that characterize the repeller can be obtained from the so called topo- 
logical pressure P(/3) [36|], which can be determined from the periodic orbits that form part of the 
repeller as the leading zero of a Ruelle (^-function 

n ('-txiw)- ^ 

p=p.p.o \ \ P\ / 

The mean Lyapunov exponent for the trajectories on the fractal repeller is A = —dP(l)/d/3, 
the escape rate is 7 = — P(l) and the topological entropy is given by ht op = P{0). In two- 



dimensional Hamiltonian systems the Kolmogorov-Sinai entropy follows from hxs = A — 7 371 ] . 
The partial Hausdorff dimension is obtained from P(dn) = and the information dimension is 
given by di = 1 — 7/A. 



V. PERIODIC ORBITS AND THE REPELLER IN THE MERIDIAN PLANE (L z = 0) 

From now on we shall focus on the analysis of the escape dynamics of photons from the black 
holes in the meridian plane. As we have just mentioned, there are some important objects that 
can be used to describe this dynamics. First we need to identify the set of periodic orbits and their 
corresponding primitive periods. To analyze the linear stability of each unstable periodic orbit we 
must integrate its fundamental matrix in time up to one primitive period. The eigenvalues of the 
resulting matrix give the stretching factors associated with the orbit. Then, once we know the 
periods and the stretching factors of the unstable periodic orbits, the main quantities of chaos that 
characterize the dynamics can be obtained following the methods described in the previous section. 

The restriction on the motion to the meridian plane L z = simplifies the set of equations to be 
integrated. Without loss of generality we can consider that the photons move in the {x, z)— plane. 
Then the equations of motion ([7]) with respect to the affine parameter [i reduce to 

a = —U~ 2 p a p a = —2Ud a U (a = x,z), (34) 

being the equations ([8]) with respect to the killing time t, 

a' = -U~ 4 Pa p' a = -2U~ l d a U {a = x, z). (35) 

These equations of motion and the equations associated with the stability analysis were simul- 
taneously integrated using a Runge-Kutta-Fehlberg method. The search for periodic orbits was 
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performed using the Newton-Raphson method. As we just mentioned, to compute the stretching 
factors of the set of unstable periodic orbits, each orbit was numerically integrated together with 
its corresponding fundamental matrix up to one primitive period of the orbit. Considering that 
the dynamics is highly unstable, it is important to compute the unstable periodic orbit with high 
accuracy in order to obtain the stability factors with good precision. In our numerical results the 
required high accuracy in the period of the orbits and the stretching factors was achieved by means 
of a Mathematica code in which a very high precision was demanded. 

The periodic orbits in the system can be classified according to a symbolic coding. In this work, 
taking as reference the symmetry axis that contains the two black holes, in our coordinate system 
the z-axis, we introduce a series of symbols to label these orbits. An orbit can cross the z-axis 
either above the black-hole located at z = Zbh, below the black-hole located at z = —z^ or in 
between the two black-holes. We label each crossing that occurs at z > z^h with the symbol "+", 
to indicate a crossing at z < —Zbh we use the symbol "— ". When an orbit crosses the z-axis twice 
at the same point in between the two black holes, drawing an "x" on the axis, we use the symbol 
"o". Any other crossing in between the two black holes is indicated with the symbol "+" or the 
symbol "— " depending on the sign of the z coordinate at the crossing point, z > (+) or z < (— ). 
Figure (JTJ) displays some of the unstable periodic orbits numerically found and their corresponding 
symbolic coding. The primitive period of each orbit is given in Table (jlj). 



In 24| a different coding to label the periodic orbits was introduced. This coding was used to 
compute the topological entropy per unit symbol. We remark that one of our goals in this work is 
to compute the topological entropy per unit time. 

In Table (|T|) are listed the non-trivial stretching factors of the periodic orbits depicted in figure 
(PQ). We recall that these are the stretching factors that are related to the directions transverse to 

the periodic orbit. The trajectories {++}({ }), {H — } and {+ o — h} have the lowest stability 

eigenvalues. The rest of the orbits are much more unstable as they have larger stretching factors. 
In general the more stable periodic orbits should dominate the asymptotic escape of particles from 
the two black-hole configuration. 

In the two black-hole system the set of unstable periodic orbits are the only critical elements of 
the flow in the fractal repeller that determines the escape dynamics of photons. This distinguished 
subset of solutions plays a crucial role as the averages of the dynamical quantities can be obtained 
in terms of periodic orbits, provide that it is dense on the repeller 211 ] . 



The dynamics in the proximity of the invariant set that defines the repeller can be analyzed 
considering the time interval that it takes to an orbit to escape from a finite spatial domain U, 
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Symbolic coding 






Period T p (/j.) 


Period T p (t) 


Stretching factor A p 


{++},{—} 






4.1180374 


37.9255026 


99.1244899 


{+-} 






15.292055 


52.7275448 


215.721472 


{+0-0} 






10.608412 


79.8770864 


16996.3144 


{+-++},{- + —} 






18.900782 


91.7239060 


30044.8517 


{+ + + 0-0}, { 04 


-0} 




14.739760 


117.837331 


1.70990814 x 10 6 


{+ + + -—} 






22.635603 


130.594585 


3.97082132 x 10 6 


{H ho -0}, {- + - 4 


-0} 




25.410381 


133.342257 


4.57097810 x 10 6 


{+- + + +-U- + -- 


-+} 




34.185532 


144.455986 


6.49738275 x 10 6 


{+ + + 0} 






18.871204 


155.798106 


1.72085241 x 10 8 


{+ ++},{- + 


+ + - 


f + — } 


26.748625 


168.528893 


3.95770189 x 10 8 


{+- + ° 0}, {- + 


- + 


+ +0} 


29.542461 


171.300173 


4.59225358 x 10 8 


{H ho +-}, { 


- + - 


+ + + 0-+} 


44.828729 


224.030938 


9.92347816 x 10 10 



TABLE I: The periods and the stretching factors of the unstable periodic orbits depicted in figure (p}. The 
periods corresponding to both the affine parameter fj, and the killing time t are indicated. 

when evolves both forwards and backwards in time. To numerically study the escape dynamics we 
introduce the forward and backward escape time functions 



211 ]. which are denned by: 



T+(;c) = Max{T >0; $^GW, V/zG[0,T[, Vie U] 
Tu(x) = Min{T < ; fteW, V// € ]T, 0] , ViG U] 



(36) 
(37) 

. The ensemble 



The escape time functions provide an useful method to construct the repeller 
of initial conditions whose escape time coordinate is larger than a given value T contains all the 
trajectories that still remain inside the domain U at time T. In the long time limit this set contains 
the trapped trajectories of the repeller and its stable manifolds in U. As time evolves this set splits 
into smaller pieces and eventually becomes a fractal set. Similarly it is possible to find the set of 
initial conditions that have already arrived in the domain U before some time —T in the past. The 
repeller is given by the intersection of these two sets. 

The compound escape function Ty(x) + |7^7(a;)| is very useful to illustrate the structure of 
the repeller. To numerically compute this function we analyze the evolution, both backwards 
and forwards in time, according to the equations of motion ([7]), of a large ensemble of photons 
initially distributed at random in a finite region of space and certain fixed initial directions. In 
order to distinguish the possible outcomes of the photons in the two black-hole field, we color 
code the initial position black if the corresponding photon falls into any of the two black holes 
and gray if the photon escapes towards infinity. Figure ([2]) displays a portrait of the basins of 



15 





A 


7 


has 


htop 


d H 


d T 


Affine parameter fi 


0.364 


0.348 


0.016 


0.162 


0.221 


0.045 


Killing time t 


0.111 


0.089 


0.022 


0.024 


0.221 


0.196 



TABLE II: Characteristic quantities of chaos related to the fractal repeller associated with the unstable 
periodic orbits given in figure (HJ. A is the mean Lyapunov exponent of the repeller, 7 the escape rate, hxs 
the Kolmogorov-Sinai entropy, h top the topological entropy, dn the Hausdorff dimension and dj the partial 
information dimension. 

attraction associated with these two possible outcomes in the plane (x, z) of initial positions, for 
an ensemble of photons with initial directions randomly fixed between five selected values for the 
angle with respect to the x-axis. It is well known that for chaotic systems the interweaving of 
outcomes leads to fractal basin boundaries; and in the two black-hole system, as expected, we can 
see them explicitly. The blown up regions depicted in the different panels of figure (J3j) exemplify 
the repeated fractal structure in the basin boundaries that characterize the hyperbolic repeller. In 
the following sections we compute different characteristic quantities of the dynamics related to this 
fractal repeller. 



VI. ESCAPE RATE IN THE MERIDIAN PLANE AND OTHER CHARACTERISTIC 

QUANTITIES OF CHAOS 



Using the primitive periods T p and stretching factors A p of the set of unstable periodic orbits 
we can compute the escape rate of the system as well as other characteristic quantities of the 
dynamics related to the fractal repeller, following the methods described in Section IV. To this 



aim we consider the cycle expansion method 
hyperbolic systems 



341 ] which has been successfully applied to other 



301 ] . Due to the high instability of higher order periodic orbits (longer symbolic 

coding) only three main periodic orbits, {++}({ }), {H — } and {+0 — |-}, contribute significantly 

to the escape rate. Indeed they give an acceptable escape rate and the different quantities of chaos. 

Figure (J4j) shows the pressure function obtained from these three periodic orbits. The main 
quantities of chaos obtained from this pressure function and its first derivative are listed in Table 
([LT1) . In particular, the escape rate which is equal to 7^ = 0.348 for the affine parameter and 
7i = 0.089 for the killing time. From the analysis of the linear stability of the set of unstable 
periodic orbits, it could be assumed that the escape dynamics of photons from the two black holes 
is mainly controlled by the outermost periodic orbit {+— }. If only this orbit is considered it follows 
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an escape rate 

r, i In I A i I (i \ In I A_i I , , 

} ~ itti = °' 351 7 ' = itw = 0102 (38) 

which is close to the escape rate obtained from the pressure function. Thus, we can conclude that 
the orbit {+— } dominates the effective escape of photons from the two black-hole configuration, 
and even provides a fair estimate of the escape rate. 

In the following section we contrast the escape rate value derived from the analysis of the 
linear stability of the periodic orbits with the escape rate obtained from numerical simulations 
that consider the time evolution of an statistical ensemble of photons. 



A. Numerical simulations with statistical ensembles 



We consider an ensemble of photons with initial positions uniformly distributed along the z-axis 
and initial velocity parallel to the x-axis. This together with the energy conservation completely 
defines the ensemble. Once all the initial conditions have been fixed the time evolution of each 
photon of the ensemble is numerically integrated. Some of the photons fall into one of the two 
black holes, while others escape to infinity. We are precisely interested in the second subset of 
trajectories. To compute the escape rate we consider a circle of control that covers the two black 
holes and which is larger than the outermost periodic orbit {+— }. Hence all the photons that cross 
out this circle escape towards infinity and do not come back to the proximity of the black holes. 

Let us denote by N(t), with r the affine parameter fi or the killing time t, the number of 
photons that remain inside the circle of control at time r (this includes the photons that fall into 
the black-holes). This function decreases with r until it reaches an asymptotic value iV(oo), which 
gives the number of photons that fall into the black holes. The escape rate 7 is defined by the 
exponential decrease to zero of the function [N(r) — iV(oo)] ~ [-^(0) — N(oo)]e~' yT as r —* 00. 

Figure © shows the results obtained from different simulations that include an increasing 
number of particles in the statistical ensemble. The larger is the number of particles in the ensemble 
the longer persists the exponential decrease with r. From the slope of these curves we extract a 
numerical escape rate 7^ = 0.344 for the affine parameter and ^ num ) — 0.083 for the killing 
time. The agreement with the escape rate obtained from the analysis of the linear stability of 
the unstable periodic orbits is excellent, see Table (jTlj). This result confirms our expectation with 
respect to the subset of periodic orbits that plays the main role in the escape dynamics. The small 
oscillatory component in the evolution of the function [N(r) — iV(oo)] can be explained by the 
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lower resonances of the system. This point has been emphasized in 30|. To further illustrate this 



( num ) 

aspect we have computed the power spectra of the numerical data [N{fi) — iV(oo)]e 7 ^ ^ for the 

(num) , 

affine parameter, and [N(t) — A^(oo)]e 7 * for the killing time. The figures (J6|) and (|7]) show a 
peak structure associated with the oscillations of both functions once the main exponential decay 
has been subtracted. The positions of the peaks is compared with the locations of the imaginary 
part of the lower resonances computed from the zeta function taking orbits that include up to six 
elements in the symbolic code. The resonance distribution explains the peak structure observed in 
the power spectrum. 

The escape rate of photons depends on the masses of the black holes and the separation distance 
between them. The larger are the masses the more extensive in space are the unstable periodic 
orbits that form part of the fractal repeller that marks the boundary between dynamical stability 
and instability, and therefore the regions where the light rays fall into the black holes. It could be 
expected then that the escape rate of photons decreases as the masses of the black holes become 
larger. Figure ([8]) illustrates this behavior and shows the decrease of the escape rate as the mass 
of one of the black holes is increased. Figure ([9]) shows the dependence of the escape rate on the 
separation distance between the two black holes. In both cases, after an initial transient, the escape 
rate decays asymptotically with the distance d. A similar behavior is observed in other scattering 
systems where the escape rate decays as f(d)/d, being f(d) a function determined by the dispersion 



of the scattering interaction potential and 1/d is the contribution of the time of flight 2l|, 1301 ] - 

In figure ([9]) we consider the scattering of photons from two black holes with masses M\ = M2 = 
M, located at different positions (0, ±d). From and (|7|8p it can be easily shown that a system of 
two black holes with scaled masses fh\ = fhi = M/d located at (0, ±1) satisfies identical dynamics 
with respect to an scaled time f defined by da/dr = da/df (df = dr/d), being a = a/d(a = x,z) 
scaled coordinates. Hence if we denote 7^ the escape rate from two identical black holes of mass M 
located at (0, ±<i), and 7 the escape rate from two black holes of the same mass fh = M/d located 
at (0, ±1), it follows that 7^ = j/d. 

VII. ESCAPE DYNAMICS OF PHOTONS FROM A PERTURBED BLACK HOLE 



In this section we study the escape rate from an extreme Reissner-Nordstrom black hole of mass 
Mi that is slightly perturbed by the interaction with an extreme Reissner-Nordstrom black hole of 
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mass M 2 located at distance d. In polar coordinates, the Hamiltonian can be written as 



n 



with the function U given by: 



h 

IP 



1 

IP 



U 



, Mi 
1 + — + 



1 / 2 , y 

2 \ Pr + r 



El 
2 



(39) 



(40) 



VP + d 2 — 2rd sin 

The analysis of the dynamics of this system can be greatly simplified if we consider the evolution 
with respect to a time v defined by the transformation 



d_ 

dfi 



1 d 
TPdv' 



(41) 



Thus we focus on the equations of motion associated with the Hamiltonian h, 

dX v dh via \ 

^ = £.— X = (r,9, Pr , Pe ) 



(42) 



In these equations, h can be analyzed as a kinetic energy term plus an interaction term — [7 4 /2. 
The function U 4 may be expressed as an expansion in powers of (d/r), with the first two terms 
given by: 



IP 



1 H H ^— 1 H sin 6* + O 



(43) 



r \ r 

In the limit (Mi + M 2 ) 2 » 2M 2 d the dynamics of the system is dominated by the radial term 
in the expansion (143 h . which is perturbed by the remaining angular terms. Thus to zero order 
approximation the space time associated with the two extreme black holes becomes spherically 
symmetric. There exists a single unstable periodic orbit which is a circle of radius r = M = 
(Mi + M 2 ) and period T v = Mtt/2. The time evolution along this orbit is given by 



r s {y) =M 

Prs{v) = 



e a {y) -- 

P6s{v) 



9 +^ 
9o+ M 

L 



(44) 



with angular momentum L = 4M. The periodic orbits associated with the affine parameter and 
the killing time are also a circle of radius r = M, but with periods Tn = 2irM and Tf = 8ttM 
respectively. 

To zero order approximation the escape rate of photons from an effective black hole of mass 
M can be determined from the analysis of the linear stability of the circular orbit (|44p . X s = 
(r s ,6 s ,p rs ,pQ s ). This implies the study of the time evolution of a small perturbation, 5X = 
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(5r,56,5p r ,5pQ), with respect X s . In terms of the Hamiltonian h, the linear differential equations 
that describe this evolution, d(5X)/dv = \(u) ■ 5X with = £ • J^r|x s , are: 





' Sr \ 




d 


56 




dv 


5p r 






V s Pe ) 





10 







M- 



-At -At 







^ <5r ^ 















(45) 



and the fundamental matrix M(X s ,i/), solution of M'(X s ,v) = \(v) ■ U(X s ,v) with M(X s ,0) = 1, 
is given by 



M(X s ,u) 



C{v) 







\C(y) - 1] 




1 


[1 - C{v)] 









C(v) 













1 



(46) 



with 



C(y) = cosh ( ^^-v 



M 



and S^z/) = sinh ( ^-r-v 



M 



(47) 



The stretching factors associated with the circular orbit are obtained from the eigenvalues of the 
matrix M(X S , v), 



1, 1, e m v , e m 



evaluated at one primitive period T v . That is, 

A={l,l, e -^, e ^}. 



(48) 



(49) 



We recall that these stretching factors are invariant under time transformations, and therefore are 
identical for the affine parameter and the killing time. As expected, two of them are equal to one. 
The values e~^ n and e v ^ w are related to the stable and unstable manifolds of the unique unstable 
periodic orbit respectively. The Lyapunov exponents associated with this orbit are given by 

lnlAI 



A,, 



1 11 



(50) 



in the time evolution dictated by the affine parameter, and 

In |A| f_ _ 1 1 



A, 



0,0, 



(51) 



20 



in the evolution with respect to the killing time. Thus to zero order the escape rate of photons 
from the unstable circular periodic orbit around the effective extreme black hole of mass M is 
determined by the leading exponent yjp = l/\/2M and y^ = 1/A\/2M = 7^/4. 

In this system there is a unique unstable periodic. Since the dynamics on the repeller is regular, 
the Kolmogorov-Sinai entropy is zero, h^s = A — = 0. Thus the system is hyperbolic but 
nonchaotic. The Pollicott-Ruelle resonances a pr are given by the zeroes of the Zeta function (|30|) 
associated with the single unstable periodic orbit, 

Z(a pr ) = ]J ( 1 ~ TaTaF ) =0, (52) 



being T the period Tn or T t . That is, 



a M(k , n ) = -it+2l + i " (53) 



for the affine parameter and 



(jW (k, n) = — _ ' + i-rv7 ( 54 ) 
p A^/2M AM v ' 

for the killing time, with k = 0, 1, 2, ... and n = 0, ±1, ±2, ... The resonances belong to the 



ower 



half plane of the complex plane a (Re a pr < 0) and, as occurs in the two-disk scatterer 2l|, |39( |. 
their spectrum forms a semi-infinity periodic array. Here the spacing along the (Re a)-axis is given 
by the escape rate 7^ 0) = l/y/2M (^ 0) = 1/4 V / 2M), and by the frequency = 1/M (w t = 1 /AM) 
along the (Im cr)-axis. These complex resonances play an important role in the time evolution of an 
statistical ensemble of photons; they determine the different decay modes and their frequencies in a 
typical scattering process. The ensemble dynamics is ruled by the resonances which are the closest 
to the imaginary axis (k = 0). The real part of these leading resonances controls the exponential 
decay on the longest time scale, which defines the escape rate of the system; and their imaginary 
parts give the frequencies of the oscillations that appear superimposed on the gross exponential 
decay. 

We now study the effect of the first angular term in the expansion (I43p on the zero-order circular 
orbit (|44j) associated with the leading radial term. We still assume the limit M 2 = (Mi + M2) 2 >> 
2M2<i, in which this term can be treated as an perturbation to the spherically symmetric motion. 
To make a first order perturbative analysis we write the U 4 function in the interaction term in the 
form 

4 a Ti/r^A / n/r\ 3 



rri (, M \ AM 2 d ( MY . / 

u =( 1+ v) + sme ( 55 ) 
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where we have introduced a perturbative parameter e, which is set equal to one at the end of the 
analysis. The perturbed trajectory can be written as 

r\{u) = r s {v) + er c (u) , p r i(v) = Prs(v) + ep rc (y) 

6i{v) = 6 s (v) + eO c (v) , po\{v) = p es {v) + epedv) (56) 



with first-order corrections to the circular orbit (|44j) that satisfy the equations 

dr c _ 

, — Pre 

dv 

de c 1 

^7 = F {p0c ~ 8rc) 

dp rc 56M 2 d {4v\ 8 

16M2d sin^V (57) 



The periodic solution of these equations gives the first-order perturbed orbit, 

r l( ,)=M + — cos(-j 
„ , , n 4 c?M 2 . (Av 

4dM 2 . [Au 



mM = 4M + ^cos(-j. (58) 

This is a quasi-circular orbit that nearly reproduces the periodic orbit {+— } in two black- 
hole systems in which {M\ + M 2 ) > d, and the periodic orbit { } ({++}) in systems where 

(Mi + M 2 ) < d, see figure (llOp . The more complex periodic orbits (longer symbolic coding) are 
out of scope of our perturbative analysis. However, in the limit (Mi + M 2 ) 2 » 2M 2 d these 
are highly unstable periodic orbits that play a secondary role in the escape rate of photons from 
the black holes. Here the escape rate is controlled by the quasi-circular periodic orbit, {H — } or 

{ } ({++}), which results from the slight deformation of the circular orbit of radius M. Hence 

the analysis of the linear stability of the zero- and first-order periodic orbits, (|44p and (|58p . should 
provide a good estimate of the escape rate of photons from a perturbed black hole. Indeed, as 
the results in Table (jlllj) indicate, both the analytical escape rate ( 7^ ) derived from the 
analysis of the linear stability of the zero-order circular orbit (|44p and the escape rate 7^ ( 7^ ) 
obtained from the leading eigenvalue of the fundamental matrix associated with the first-order 
quasi-circular orbit (|58l) are very close to the escape rate 7^ niim ^ ( ^ num ^ ) calculated from the 
numerical simulations that consider an statistical ensemble of photons, see figure (jlip . As figure 
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System 


(0) 




(num) 


(0) 

It 


7* (1) 


(num) 

it 


(Mi =1, M 2 = 0.01, d = 2) 


0.70011 


0.69559 


0.69258 


0.17503 


0.17390 


0.17624 


(Mi = 7, M 2 = 0.1, d = 2) 


0.09959 


0.09957 


0.09901 


0.02490 


0.02489 


0.02436 



TABLE III: Escape rate values for two systems of perturbed black holes, (7} ) is the escape rate 

derived from the zero-order circular orbit (|44p . 7^7 (7^ ) is the escape rate obtained from the analysis 
of linear stability of the first-order orbit (|5"5|) and 7^™""^ ( ^ ntlm ) ) [ s the escape rate calculated from the 
numerical simulations that consider an statistical ensemble of photons. 

([8]) shows, the zero-order analytic estimate ( 7}°' ) even provides an approximated value for 
the escape rate from a system of two black holes with similar masses. 

In the numerical simulations with an statistical ensemble of photons, see figure (jllh . the escape 
from the perturbed black hole presents a smooth oscillation superposed on the gross exponential 
decay determined by the escape rate 7^ n " (7^™""^ ). The zero-order Pollicott-Ruelle resonances 
(|53h and (I54j) account for this behavior of the decay curves; the exponential decay determined by 
the escape rate 7^ ~ 7^ ( 72 ~ 7^ ) is given by the real part of the leading resonances k = 0, 
and the frequency of the oscillation ~ 1/M {wt — 1/4M) is given by the imaginary part of 
the leading resonance given by n = 1 (A; = 0). This is illustrated in figure (|12h which displays 
the Fourier spectra of the signals in figure (jlip . It clearly shows that the most intense peaks in 
the spectra, associated with the frequency of the main oscillation in the escape rate, are localized 
pretty close to these frequencies predicted by the leading zero-order Pollicott-Ruelle resonances. 

These peaks correspond to the frequency of the leading orbits, the orbit { }({++}) i n case A 

and the orbit {-| — } in case B, see figure f j 1 [) . This is another indication of the leading role of the 
quasi-circular orbit that results from the slight deformation of the circular orbit of radius r = M 
in the escape dynamics of photons. 



VIII. CONCLUSIONS 



We have studied the scattering of photons in the Majumdar-Papapetrou static spacetime of 
two extreme Reissner-Nordstrom black holes held fixed in space due to the balance between the 
gravitational attraction of their masses and the electrostatic repulsion of their charges. We have 
identified the set of unstable periodic orbits that form part of the fractal repeller which fully 
characterizes the chaotic escape dynamics of the photons from the two black holes. These orbits 
were classified according to a symbolic coding. 
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The linear stability of the dynamics, in particular the unstable periodic orbits, was analyzed 
using their stretching factors, which are given by the eigenvalues of the fundamental matrix inte- 
grated up to one primitive period of the orbits. With the primitive periods and stretching factors 
of the periodic orbits we determined the topological pressure, and from this function and its first 
derivative the main quantities of chaos that characterize the fractal repeller. In particular the 
escape rate, which was calculated using the trace formula for hyperbolic flows of Cvitanovic and 
Eckhardt. This escape rate derived from the periodic orbit theory was in good agreement with 
the value obtained from numerical simulations that analyze the escape dynamics of an statistical 
ensemble of photons. 

From the study of the linear stability of the dynamics we also identified the periodic orbit that 
plays the leading role in the escape of photons from the two black-hole configuration. In systems 
with two identical extreme Reissner-Nordstrom black holes the escape is mainly controlled by the 
outermost periodic orbit {-| — } that encircles the two black holes. In systems where a main black 
hole of mass Mi is perturbed by the weak interaction with a black hole of mass M2, the escape 
of photons is dominated by the quasi-circular orbit that results from the slight deformation of the 
circular orbit of radius r = M = (Mi + M2). Depending on the separation distance between the 
two black holes this leading orbit can cover the whole system, and corresponds to the periodic orbit 
{-I — }, or only encircles the main black hole, and is given by the periodic orbit {+- 1-} ({ })• 

In systems of a perturbed black hole, the analysis of the linear stability of the unique zero- 
order periodic orbit in a perturbative treatment provides an analytic estimate of the escape rate, 
7^ = \j\[2M (jt^ = l/4\/2M), which is in good agreement with numerical value derived from 
the analysis of the linear stability of the first-order perturbative periodic orbit, and also with the 
escape rate obtained from the time evolution of an statistical ensemble of photons. This analytical 
value even provides an approximated value for the escape rate of photons from a system of two 
extreme Reissner-Nordstrom black holes with similar masses. And therefore a way to estimate the 
mass of a black hole from the escape rate of photons. 
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{+0-0} {+-+0-0} {+-+0 — o} 




{+++0 — o} {+++0-0} {+-+0 — 0+-} 



FIG. 1: Some of the unstable periodic orbits in a system with two black holes located at z^h = ±1 and 
masses Mi = M 2 = 1. The symbolic coding that labels each orbit is indicated. Due to the symmetry of 
the system with respect to equatorial line (the x-axis), the change of sign of the coordinate z in each orbit 
on the last two columns gives an additional periodic orbit. The symbolic coding for these inverse orbits is 
obtained exchanging the symbols "+" and "-" in the orbits depicted in these columns. For instance, the 

inverse orbit of {++} has symbolic coding { }, the inverse of {H h+} is { — I }, the inverse of 

{H ho — o} is { — I — o +o} and so on. 
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FIG. 2: Outcomes for photons in the two-black-hole field. The black dots are the initial positions of the 
trajectories which fall into any of the two black holes. The gray dots localize the initial positions of the 
photons that escape towards infinity. The evolution according to the equations of motion ([7]) was considered, 
both backwards and forwards in time. The initial positions (x, z) were chosen at random in the finite region 
depicted, and the initial directions were randomly selected with an angle with respect to the x-axis equals 
to 0, tt/5, 27r/5, 3tt/5 or 4tt/5. 
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FIG. 3: Several blown up regions of figure ||3J) 



showing the fractal structure of the repeller. 
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FIG. 4: The pressure function (solid line) obtained considering the periodic orbits that include up to six 
elements in the symbolic code, see figure ([T]). The dashed line is the first derivative of this pressure function. 
The panel (a) is the pressure associated with the affine parameter, and (b) the pressure corresponding to 
the killing time . 




FIG. 5: The evolution of [N(n) — N(oo)] ( [N(t) — N(oo)) ) with the affine parameter (killing time) for three 
ensembles with an increasing number of photons, 5 x 10 5 , 5 x 10 6 and 5 x 10 7 . The dashed lines give the 
best linear fitting for the ensemble with the highest number of photons. 
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FIG. 6: In the upper panel the oscillatory component of the function [N(fi) — N(oo)} depicted in the left 
panel of figure (O, once the main exponential decay has been subtracted. In the middle panel its power 
spectrum, being the vertical dotted lines the imaginary part of the resonances given in the lower panel. 




FIG. 7: The same as figure (|6]), but for the function [N(t] — iV(oo)] ) depicted in the right panel of figure 
©■ 
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FIG. 8: The escape rate 7^ ( j t ) from the two black holes located at Zbh — ±1 for different values of the mass 
Mi (M 2 = 1)- The data were obtained considering ensembles of photons with initial velocity parallel to the 
x-axis and initial positions uniformly distributed along an interval on the z-axis. The circles correspond to 
an ensemble initially distributed in the interval [1-1,4] and the crosses to an ensemble initially distributed 
in the interval [—4, —1.1]. The dashed lines are the analytic estimate of the escape rate from a perturbed 
black hole, 7 ^ 0) = 1/V^(M 1 + M 2 ) and 7t (0) = l/A^{M 1 + M 2 ), see Section VII. 
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FIG. 9: Escape rate 7^ ( 7* ) from two identical black holes (Mi = M 2 = 1) for different separation distance 
between them. The black holes are located at Zbh = ±<^- The data were obtained using ensembles of photons 
with initial velocity parallel to the x-axis and initial positions uniformly distributed along the z-axis. 
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FIG. 10: Some of the periodic orbits in two systems of a black hole of mass M\ perturbed by the interaction 

with a black hole of mass M 2 located at distance d. The solid lines give the periodic orbits {++}, { }, 

{H — }, {H h+} and {+ o — o}, and the dashed line is the first-order periodic orbit (|58|) . The system (a) 

corresponds to (Mi = 1, M 2 = 0.01, d = 2) and (b) to (Mi = 7, M 2 = 0.1, d = 2). In both systems the 
main black hole of mass Mi is located at (0, 0) and the black hole of mass M 2 at (0, d). 
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FIG. 11: The evolution of [N(fi) — iV(oo)] ([N(t] — JV(oo)]) with the affine parameter (killing time) for two 
statistical ensembles in the perturbed systems of figure (jTTJJ) . Panels (a) and (c) correspond to the system 
(a) in figure (TIT)]) , and panels (b) and (d) to the system (6) in figure (flU)) . The solid line gives the escape 
from a circle of control centered in between the two black holes, at (0, d/2), and the dashed lines the escape 
from a circle of control centered on the main black hole, at (0,0). See the Section VII for comments. 




FIG. 12: The Fourier transform of the four signals in figure (fTTj) . The dashed lines localize the frequency 
of the main oscillation in the escape rate predicted by the leading zero-order Pollicott-Ruelle resonances, 
UV =s l/M (w t ~ I /AM). 



